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ABSTRACT 

I consider the nonaxisynimetric linear theory of a rotating, isothermal magnetohydrodynamic 
(MHD) shear flow. The analysis is performed in the shearing box, a local model of a thin disk, using a 
decomposition in terms of shearing waves, i.e., plane waves in a frame comoving with the shear. These 
waves do not have a definite frequency as in a normal mode decomposition, and numerical integra- 
tion of a coupled set of amplitude equations is required to characterize their time dependence. Their 
generic time dependence, however, is oscillatory with slowly-varying frequency and amplitude, and 
one can construct accurate analytical solutions by applying the Wentzel-Kramers-Brillouin method 
to the full set of amplitude equations. The solutions have the following properties: 1) Their accuracy 
increases with wavenumber, so that most perturbations that fit within the disk are well-approximated 
as modes with time-dependent frequencies and amplitudes. 2) They can be broadly classed as in- 
compressive and compressive perturbations, the former including the nonaxisymmetric extension of 
magnetorotationally unstable modes, and the latter being the extension of fast and slow modes to a 
differentially-rotating medium. 3) Wave action is conserved, implying that their energy varies with 
frequency. 4) Their shear stress is proportional to the slope of their frequency, so that they transport 
angular momentum outward (inward) when their frequency increases (decreases) . The complete set of 
solutions constitutes a comprehensive linear test suite for numerical MHD algorithms that incorporate 
a background shear fiow. I conclude with a brief discussion of possible astrophysical applications. 
Subject headings: accretion, accretion disks — MHD 



1. INTRODUCTION 

The combined effects of magnetic fields and shear are 
important i n such astrophysical systems as gal actic disks 
(IShul 11974 . accretion disks (iBialbus fc Haw lcv 199li 
11998 ) and differentiallv-rotating stars (jSpruiti i200^ . 
A simple theoretical model for understanding the lo- 
cal behavior of a magnetohydrodynamic (MHD) shear 
flow is the shearing box, developed originally for the 
study of gravitational spiral \ yaves in galactic disks 
(jGoldreich fc Lvnden-Bellfl965l : rhilian fc Toomrelll966l ^ 
and used extensively in accretion disk theory. The model, 
which is local but fully nonlinear, is generated by an ex- 
pansion of the equations of motion in the ratio of the 
disk scale height to the disk radius, and is therefore rig- 
orously valid only for a thin disk. A great deal of in- 
sight into the behavior of astrophysical disks has been 
obtained through both analytical and numerical studies 
of the shearing box. This paper extends previous work 
on the linear theory of the MHD shearing box, with the 
primary goal of obtaining analytical solutions for use in 
testing numerical algorithms. 

A linear analysis in the shearing box is complicated 
by the presence of the shear, the effect of which upon a 
plane wave of definite wavenumber is to rotate its wave 
crests into a trailing configuration with a primarily ra- 
dial wavevector. The t ime dependence of these shear- 
ing waves, or "shwaves" (! Johnson fc Gammiell2005l ) . can- 
not in general be expressed in terms of a frequency 
eigenvalue as in a normal mode decomposition; one 
must solve instead a coupled set of differential equa- 
tions. As a result, most linear analyses in the shear- 
ing box consider only axisymmetric modes, the spa- 
tial configurations of which remain stationary in the 
presence of shear. Characterization of the nonaxisym- 



metric shwaves generally requires numerical integration 
of the set of equations governing their time depen- 
dence (e.g.. [Bal bus & Hawlcv 1992; Chagclishvili et al.l 
Il997t iFan fc Loulll997t iBrandenburg fc DintransI 12006V 
although approximate soluti ons have been o b tained 
in some limiting cases. IBalbus fc HawievI (Il99l 
demonstrate the transient amplification of incompressive 
shwaves in the presence of a destabilizing field and show 
that extending the axisymmetric dispersion relation to 
general wavenumber provides a good approximation to 
the growth rate d uring the phase of ne arly exponen- 
tial amplification. 1 iKim fc Ostriked ()2000l ) analyze MHD 
shwaves in the limit of small shear to obtain what they 
refer to as "coherent wavelet solutions" ; these modal so- 
lutions are valid over a time scale that is short compared 
to the shearing time scale. 

A glance at the numerical evolution of a set of 
shwaves, however, clearly indicates that their generic 
time dependence is oscillatory with a slowly varying fre- 
quency and amplitude, even in the presence of significant 
shear. Such a time dependence ought to be amenable 
to a Wen tzel-Kramers-Brillouin ( WK B) analysis; in- 
deed^^both IBalbus fc HawlevI (|1992f ) and lKim fc Ostriked 
(|2000f ) note that their dispersion relation represents the 
lowest order solution in a WKB expansion. The purpose 
of this paper is to generalize these analytical results into 
uniformly valid asymptotic solutions for both compres- 
sive and incompressive shwaves and for arbitrary values 
of the shear, and in particular to take the expansions to 
the next order and solve for the slowly-varying ampli- 
tudes. 

The WKB method is typically applied to a second- 



^ A similar analysis is performed bv lRvu fc GoodmanI ((1992) in 
the context of convective instability in disks. 
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order differential equation, but the method can be ap- 
plied in a straightforward manner to a coupled set of dif - 
ferential equations of arbitrary order (|Brethertonlll968l ). 
With time as the independent variable, the WKB phase 
can be regarded as an integral over a local ("instanta- 
neous") frequency. The lowest order terms thus yield a 
time-dependent dispersion relation, a solution of which 
can be integrated to obtain the WKB phase (the eikonal 
equation) . The next order terms yield a first-order differ- 
ential equation for the slowly-varying amplitude of one 
of the dependent variables (the transport equation) . The 
amplitudes of the other variables are then determined by 
the eigenvector components from the lowest order ex- 
pansion. The problem is thus reduced to a set of decou- 
pled equations for a slowly varying amplitude and phase 
which can be integrated by quadrature.^ Alternatively, 
the overall amplitude can be determined from the con- 
straint that wave action is conserved, a general result for 
oscillatory solutions derived from a Lagrangian. 

I begin in fj2]with an overview of the local equations for 
an isothermal MHD shear flow in their general and lin- 
earized forms, as well as the decomposition of the latter 
in terms of shwaves. [J3]and ^describe the leading order 
solutions for incompressive and compressive shwaves, re- 
spectively, and ^J5] discusses their energy and angular mo- 
mentum transport properties. I summarize and discuss 
applications in ^J6l 

2. BASIC EQUATIONS 

As discussed in the Introduction, the shearing box is 
generated by an expansion of the equations of motion in 
the ratio of the disk scale height to the disk radius. The 
resulting model is a local Cartesian frame in which the 
orbital velocity varies linearly with radius. The shearing 
box equations in isothermal MHD are 



| + v.(H = o, 



(1) 



at p 



B VB 



2nxv- 2q9?xx = 0, 



local Coriolis force. The coordinates x^y and z corre- 
spond to the radial, azimuthal and vertical directions, 
respectively. 

The equilibrium state is one of constant density pQ, 
spatially constant magnetic field B^ and orbital velocity 



■Wo 



-q^lx y. 



(6) 



In linearizing about this equilibrium state, it is conve- 
nient to work in units in which the mean density po is 
unity, so that the magnetic field has units of velocity. In 
these units. 



p = 1 + ^p, 



B 



VA + SB, 



(7) 
(8) 



(9) 



where va = -Bo/v47r is the Alfven velocity. 

In the absence of fixed boundaries, the natural basis 
for local perturbations in the shearing box is a decom- 
position in terms of plane waves in shearing coordinates, 
{x',y',z') EE {x,y ~ V()t,z): 

S{t) exp (ik ■ x) — d{t) exp (iA:a;[0]a; -I- iky[y — vot] + ikzz) . 

(10) 

A non-shearing observer sees these as plane waves with 
a time-dependent radial wavenumber: 



kx{t) = kx{0) + qflkyt. 



(11) 



As a result of the time dependence of their radial gra- 
dients, nonaxisymmetric shwaves do not have a definite 
frequency as in a normal mode decomposition. One can 
readily verify from the induction equation ^ that 



VAyit) = VAy{0) - q^VAxt, 



(12) 



(2) 



so that Va also varies linearly with time in the presence 
of a radial field. 

Expanding equations to linear order in the per- 

turbation variables defined in equations ((7))- (|10p results 



and 



dB 

— - V X (i; X B) = 0, 



V • B = 0, 



(3) 



(4) 



where p is the volume density, v is the bulk fluid velocity, 
Cg is the isothermal sound speed, B is the magnetic field, 
fl = ^Iz is the rotation vector and 



1 dinn^ 

2 dlnr 



(5) 



is the shear parameter. The last term in equation ^ 
is the radial tidal force and the penultimate term is the 



^ IFoelizzo &: Taegeil 119941 . 119951 ) perform an analysis similar 
to what is described here, although their WKB expansion is per- 
formed in frequency space and they provide analytical expressions 
only for strongly leading and trailing shwaves. 



5p + ik ■ Sv ~ 0, 



Sv + ic^Spk + i{vA ■ SB)k — i{vA ■ k)SB 



(13) 



-2nSvyX + {2- q)n6vxy = 0, (14) 

SB — i{vA ■ k)Sv + i{k ■ 5v)va + qftSvAxV ~ 0, (15) 
and 

k-SB = 0, (16) 

where overdots denote time derivatives and where the 
common factor exp {ik ■ x) has been removed. Equations 
1131)- (HH) are solved in various limits in the following sec- 
tions. 

^ There is an additional vertical tidal term C^zz which I have 
ignored for simplicity; this is equivalent to considering scales much 
smaller than a disk scale height. 
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3. INCOMPRESSIVE SHWAVES 

3.1. Equations 

By definition, incompressive shwaves satisfy 

k-Sv^O (17) 

throughout their evolution. The appropriate limit for ex- 
tracting them is therefore the Boussinesq approximation, 
in which dt <C c^fc. This implies both a subthermal field 
(/3 = cl/v\ 1) and wavelengths that are short com- 
pared to the disk scale height (csfc/rj » 1). The behavior 
of these transverse shwaves can be isolated by taking the 
cross product of equation (|14p with k twice.'* Combin- 
ing the resulting equation with the induction equation 
and the divergence-free constraints on the velocity and 
magnetic field perturbations, one obtains the following 
reduced set of linear equations: 

St] — i{vA ■ k)k'^5vAx + 2nkz5LJx = 0, (18) 



1{VA ■ k)SJx - {2 - q)nh^ ^ 0, 



Sti 

6vAx - i{vA • 0: 



where 



SJx - i{vA ■ k)6ujj: - qflkzSvAx = 0, 



6LOx-i(vA-k)SJx-(2-q)nk,^^0, (19) 

(20) 
(21) 

(22) 

is proportional in two-dimensions to the z-component of 
the perturbed vorticity, 

Slu^ = kydvz — k^Svy (23) 

is the x-component of the perturbed vorticity, and 

SJ.^ = ky6vAz - kzSvAy (24) 

is proportional to the x-component of the perturbed cur- 
rent. The self-consistent density perturbation is obtained 
by taking the dot product of equation with k: 



5p- 



VA- SB .2Vt f tj; SVy ,-lky 5V.J; 



(25) 



In deriving the above equations, it is important to notice 
that 



— (k ■ Sv) • 5v 



(26) 



k ■ 5v — — (k ■ 5v) — ■ 5v w —qVlkySv^ 

due to the time dependence of k. 

3.2. Solution 

Without loss of generality, one can express each per- 
turbation variable in a WKB form: 



m 



oc 

E 

n=0 



5(")(t)exp (^-i J uj{t)dt^ , (27) 



where the WKB phase has been expressed in terms of 
an integral over a slowly varying frequency. For the in- 
compressive shwaves, it is not obvious a priori what small 
parameter to use in the expansion, but clearly a valid ex- 
pansion requires <c S^'^\ as well as dtilnuj) ^ uj 
(the adiabatic approximation).^ 

* This was pointed out to me by Yoram Lithwick. 

^ An additional requirement is that dtlnS^^^ <C i^, but since 
5'"^ is a function of uj, a slowly-varying frequency implies a slowly- 
varying amplitude. 



Expressing each perturbation variable in the form (|27p . 
the lowest-order terms from equations (fT8|) - ([2T1) result 
in a dispersion relation with th e same form as the ax- 
isymm etric dispersion relation ((Balbus &i Hawlevlll99lL 
[1993): 



fc2 



[VA ■ kf = 0, 



(28) 



where = 2(2 — q)VP' is the epicyclic frequency and 



UJ 



{vA ■ kf 



(29) 



The differences here are that k is the full three- 
dimensional wavenumber so that uj depends upon time 
via kx{t).^ The eigenvector components can be expressed 
as 



5v 



va ■ k 



■ k X 



[d^k^ 

k X X + i—TT-. X 

2Qk,Lu 



SvAx 



and 



5B 



k X \ k X X 



2^k^LL! 



5vax 

i.2 I 1,2 ■ 



(30) 



(31) 



As demonstrated explicitly in Appendix]^ the next or- 
der in the WKB expansion yields the time dependence of 
the slowly-varying amplitude. An alternative approach is 
to compute the amplitude from the conservation of wave 
action. It is well known that slowly-varying wavetrains 
in non-dissipative, continuous systems obey the follow- 
i ng conservation law to all orders in a WKB expansion 
(IWhithamI 119651 : iBretherton fc Garrett! 119681: IWhithamI 



d_ f£_ 

dt \Zi 



V 



£ 



0, 



(32) 



where £ is a phase-averaged perturbation energy, oj is 
the frequency as measured by an observer moving with 
the local velocity of the medium, Vg = 'Vj^uj is the group 
velocity, and £ /to is referred to as the wave action density. 
Since the phase-averaged energy of a shwave is spatially 
constant, expression ([32]) in this context reduces to 



d_ /£_ 

dt \Z} 



0. 



(33) 



The perturbation energy appearing in the adiabatic 
invariant £ /uj is 

£^u;^-C, (34) 

where /! is a phase- averaged perturbation Lagrangian.^ 
The energy as defined in expression (|34p is derived in 
Appendix [Cl in the incompressible limit, expression (jC6|l 
for the average shwave energy reduces to 



£^i{Lo'\^f~2qn'\U^ + \SvAf 



(35) 



where ji^p = ^ • ^* is the magnitude of the Lagrangian 
fiuid displacement ^ and |<5waP = SB -SB* . The relation 

^ Notice that while and fc both vary linearly with time, t)^ • fc 
is a constant. 

^ This is generally not equivalent to an expansion of the energy 
expressed in terms of Eulerian perturbations. 
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between ^ and the Eulerian velocity perturbation in a 
differentially-rotating medium is given by 



6v = ^ + qfl^xV- 



(36) 



Computing £ using the eigenvector relations (I30p and 
(Plj) subject to the constraint £ (x uj from ([55]) gives* 



SvAx oc 



(37) 



so that to within a multiplicative constant the amplitude 
of the incompressive shwaves is given by 



5v ^ uj 



and 



a)2fc2 

fc X ( fc X a; + «— x] (38) 



oB = — a; 



fex I fe X a; + I ^ — X 

(39) 

Notice that these are nonlinear solutions to the basic 
equations, a general res ult for a single incom pressive 
plane wave perturbation ([Goodman fc Xulll994f ). 

3.3. Validity of Expansion 
For wavenumbers with va • fc S> f^, both branches of 



the dispersion relation (|28p reduce to 



(40) 

Notice that these short-wavelength modes do not have 
fc II t;^ (i-e., they are not shear Alfvenic), since 6{t) 
Tr/2 due to the shear, where 9 is the angle between va 
and fe. 

For wavenumbers with va ■ k an accurate expan- 



3> ky. One way to see this is 
II21I1 in terms of Svx, which in- 



sion generally requires /c^ 
to rewrite equations (fTSll - 
troduces an additional term (2q^kxky)Svx into equation 
(fT8|) . The resulting dispersion relation matches ([28]) in 
the limit kz 3> ky (for q ~ 1). 
For wavenumbers satisfying 

len^k^ (va ■ k]^ 



< 1, 



the two branches of the dispersion relation 
to 



and 



u;l = {va -kYil 



417^ 



(41) 
reduce 
(42) 

(43) 



Notice that lo^ < for q < 0, reflecting the un- 
stable na ture of this br anch of the dispersion relation 
(Balbus fc HawlevI I1992D . The vahdity of expression 
( 4T|) generally requires va ■ k ^ 51; this implies either 
Va • fc « or, since the wavenumbers must also satisfy 
the Boussinesq approximation, ^{va/cs) <C VAk <C fl, 
i.e., /3^/2 ^ Csk/fl ^ 1. In the limit of a very weak 
field or for wavevectors oriented nearly perpendicular to 
Va, there is thus a range of wavenumbers satisfying (j4ip 



This result also requires taking the limit kz ky; see 



that increases with decreasing field strength. Eventually 
(|4T|) breaks down, however, due to the increase of k with 
time, and the frequency will asymptote to va • k. 
Equation (|42p satisfies the adiabatic approximation for 



q^ky \ kx 



< 1. 



(44) 



This ratio is a (constant) maximum at large |A:2:(t)|, and 
therefore a uniformly valid expansion requires 



> — ■ 



(45) 



These solutions have been verified with a numerical 
integration of equations ([Tn)) -(|16 p . with excellent agree- 
ment in general for wavenumbers satisfying k^ 3> ky. 
Figure [1] shows the error in the solution corresponding to 
the positive branch of the dispersion relation (|28p for 
a series of numerical integrations with varying k^/ky, 



kx{0) 



'lOky and 9 = 3/2 (Keplerian rotation). The 



vertical wavenumber is fixed at CsA:^/(27rf2) ~ 10 and the 
mean field is vaz ~ 2fl/kz. The error is defined to be 



I £n £a I 



£a 



(46) 



where £„ and £a are the numerical and analytical values 
of expression psp . A numerical check of the accuracy 
of the solutions in the limit ky 3> kz indicates that the 
errors remain small (< 30%). 

Figure [5] shows the error in the solution correspond- 
ing to the negative branch of the dispersion relation 
(|28p for a series of numerical integrations with vary- 
ing kz/ky and kx{0) = ^ky. The azimuthal wavenum- 
ber is fixed at Csky/(2Trn) = 1 and the mean field is 

VAz = \/'^5/16il/kz, corresponding to the most unsta- 
ble wavenumber in axisymmetry. 
Notice that for wavenumbers satisfying 



1,2 

Ky 



(vA ■ kf 



2q9? - {va ■ ky 



(47) 



the lower branch of the dispersion relation (|28|) goes to 
zero at 



kx {i) -^ky 



\ 



kl / 2qQ.^ 



V V>A ■ kf 



-1-1. 



(48) 



While the solution is accurate in both the oscillatory and 
exponential regimes, it breaks down as w — *■ 0, and a 
uniformly valid expansion would require a more careful 
treatment of the solution near these turning points. 

A numerical check of the incompressive solutions at 
varying values of Cgk/Q indicates that compressibility ef- 
fects are generally negligible for Csk/{2Tr^l) > 3. 



3.4. Weak Field Limit 
For Va ■ k <^ H. and k^/ky < qil/K, condition ([42 



IS 



violated and one must revisit the Boussinesq equations 
(fT8 |) -([2T |l in that limit. For general kz/ky and va ■ k <^ 
Q, these equations can be combined into a second-order 
differential equation. 



k SvAx 



2 q kj^ 



Svax + K^klSvAx = 0, 



(49) 



which can be solved in terms of hypergeometric func- 
tions: 



5vAx = CiF 



1- a 1+a 1 



4 '2' fc2 + fc2 



+C2 F 



3 — a 3 + a 3 



2' kl + k-'. 



where 



q^n^ki 



(50) 



(51) 



and ci and C2 are integration constants that depend on 
the initial conditions. The other perturbation variables 
are given by^ 



5Jx — 2Qkz J SvAx dt, 

e - 



5r] 



va ■ k 



SvAx, 



and 

with 

{6vy,6vAy) 
and 

{dv;,,SvAz) 



duJx = -« ; ovAx, 

VA ■ k 

_ kz {5lox, SJx) + kxky {dvx,SvAx) 

~ i.2 I jU2 

ky {SuJx, SJx) - kxkz {dvx,SvAx) 

1,2 I 1,2 



(52) 
(53) 

(54) 
(55) 
(56) 



3.5. Hydrodynamic Limit 

In the absence of magnetic fields, the amplitude of the 
incompressive shwaves is given by 



Sv = \/uj k X [k X X + 



and 



Csk 



l-q 



ky + k1 

fc2 



(57) 



(58) 



with Lo = nkz/k. 

These are the WKB solutions analyz e d by 
lAfshordi et all (pOOSh and iBalbus fc Hiwdevl (pOM l 
in the context of enhanced angular momentum trans- 
port in unmagnetized disks. They are uniformly valid 
for all wavenumbers satisfying condition (|45p . For gen- 
eral wavenumber, the solutions can a gain be expressed 
in ter ms of hypergeometric functions (|Balbus fc HawlevI 
[2006h . 

® The hypergeometric functions have simple integration and dif- 
ferentiation relations (e.g., Abraniowitz & Stcgun 1972) and can 
be implemented for small argument by the GNU Scientific Library. 
A Mathematica script that generates the solutions described here 
is available upon request. 

Notice that th e energy amplification of these solutions 
IIAfshordi et al.l|2005l ) is due to the time dependence of w and the 
conservation of wave action; see ^ 



4. COMPRESSIVE SHWAVES 

4.1. Equations 

In the absence of rotation, the compressive modes are 
the fast and slow waves. They have dt ^ Csk and a 
nonzero velocity component along their propagation di- 
rection. It is therefore natural to work with the linear 
equations in the following form: 

k ■ 5v + ik^ {c^Sp + Va ■ Sva) 
+ (2 - q)nkySv^ ~ 2nkjvy = 0, (59) 

Va - Sv + i {va ■ k) c^Sp 
+ {2 — q)flvAy5vx — 2Q,VAxSvy = 0, (60) 

and 

Va ■ Sva — i {va ■ k) Va ■ Sv + iv\k ■ Sv 

+qVtVAySvAx = 0, (61) 

along with the continuity equation (jl3p . 

4.2. Solution 

Establishing a uniformly valid expansion for compres- 
sive shwaves requires taking the limit Hk 3> 1, where 
H = ^Jc^^+v\/Vt\ these short wavelengths are dynami- 
cally unaffected by the rotation. The lowest-order terms 
in the expansion yield the dispersion relation 

oj^ - {cl + v\) k^uj^ + clk^ {va ■ kf = 0, (62) 
with eigenvector components 



5v = ^fc -VA-kvA]Sp 



and 



SB 



VA 



VA ■ k 



k Sp. 



(63) 



(64) 



For short-wavelength compressive shwaves, expression 
C6[) for the average shwave energy reduces to 



£^\{\Sv\' + cl\Sp\' + \SvA\' 



(65) 



where \Sv\^ = Sv ■ Sv* {— up' |^|^ to leading order in 
[_fffc]^^) and \Sp\^ = SpSp* . Computing £ using the 
eigenvector relations ([63]) and ((64|) . subject to the con- 
straint £ oc tj, gives 



Sp oc Lok 



(66) 



— c^k^ {va ■ k) 

As demonstrated explicitly in Appendix [B] the same re- 
sult is obtained by going to the next order in the WKB 
expansion. To within a multiplicative constant, then, the 
amplitude of the compressive shwaves is given by 



Sv 



SB 



ujk 



UJ 



k — Va ■ k Va 



*-{vA-kfc^^k^ \k^ 
uj'^k I u) / Va ■ k 



UJ 



UJ^ 



{va ■ k)^ c^k^ 



Va 



fc2 



and 



Sp- 



Luk 



j4 - (-y^ . fc)" C2fc2 



(67) 

(68) 
(69) 
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4.3. Validity of Expansion 



In the limit 



{vA-k)Uf3 



« 1, 



(70) 



which is vahd unless /? ~ 1, 
and slow shwaves reduce to 



the frequencies of the fast 



-/3) 



(magnetosonic shwaves) and 

c: {vA 



„2 ... . . 



ci + vi 



1 + 13 



{va ■ ky 



(71) 



(72) 



respectively. 

For a mean field that is constant with time {vax — 0), 
the fast shwaves satisfy the adiabatic approximation for 



at(lnw+) 



H 



< 1. 



(73) 



2qky 



mViikl + fc2) 



VAk ^ Cgfc and 



For /3 1, the slow shwaves have dt 
the basic assumption for compressive shwaves is no longer 
valid. The solution for the slow shwaves is therefore ac- 
curate only for /? ^ 1. Even in this regime the solution 
breaks down asymptotically in the presence of a time- 
dependent magnetic field, since in that case cx ^/P 
decreases with time; the breakdown occurs when 



UJ- 

Cgk 



VAk 



cos 6* < 1, 



(75) 



i.e., when va and k are nearly perpendicular. 

These solutions have been verified with a numerical 
integration of equations p^ - p^ . with excellent agree- 
ment in general for wavenumbers satisfying Hk /{2t:) 3> 
1. Figure [3] shows the error in the energy of a fast shwave 
for a series of numerical integrations with varying kz — ky 
and q = 3/2. The initial mean field is vax = —Cs, 
VAy = 3cs and vaz = Cs- A consistent comparison of 
the errors in this case is difficult since the frequency of 
the compressive shwaves increases with k and therefore 
the requirements for an accurate numerical integration 
become more stringent. To make a fair comparison, the 
integration time is scaled with kz and ky to give approx- 
imately the same number of oscillations during each in- 
tegration. For kz = ky, the scale factor / in Figure [3] is 
given by 



1/2 



/ 



'1 + 99 



2T:n 



(76) 



Figure m shows the error in the energy of a slow shwave 
for a series of numerical integrations with varying kz — ky 
and q = 3/2, where again the integration time has been 
scaled by /. The mean field in this case is vax = 0, 
VAy = Cs and vaz = O.lcj. 

There appear to be two additional regimes in which 
these solutions break down. Notice from equations (|63p 



and ([64)) that the compressive shwaves are singular at 
(D = , which occurs when {va ■ k) — v'^k^ cx sin^ 9 — 0, 
i.e., when k and va are parallel. There is a small subset 
of wavenumbers, therefore, for which the solution breaks 
down as the angle between k and va becomes small. 
With VAx — 0, for example, this occurs at kx{t) = 
and VAz/vAy ~ kz/ky. In addition, for wavenumbers 
satisfying expressions ([17|) and the presence of an 

unstable incompressive shwave can significantly alter the 
compress ive solutions. There is al ways coupling between 
shwaves (jGogoberidze et aDl2004[ ) since their identity as 
distinct modes is only accurate to within the errors of the 
solutions; coupling to an unstable shwave will eventually 
dominate the evolution of an initially stable shwave. 

4.4. Hydrodynamic Limit 

In the absence of magnetic fields, the amplitude of the 
compressive shwaves is given by 



This has a maximum value at k^ = {ky + k1)/2, so that and 
a uniformly valid expansion requires 



5v = \fu} k 



5p. 



(77) 



(78) 



where oj = Cgk. In the limit kz — *■ 0, these reduce 
to the WKB solution s for short-wavelength com pressive 
shwaves discussed bv lJohnson fc Gammid ()2005() . 

5. SHWAVE ENERGY AND ANGULAR MOMENTUM 
TRANSPORT 

The conservation of wave action is directly related to 
the evolution of wave energy and angular momentum 
|Balbus.2003i) . It is straightforward to show that 



Ot LO Ot 



(79) 



for both the incompressive and compressive shwaves, 
where Txy is the phase- averaged shear stress SvxSvy — 
SvAx^VAy Two things are apparent from these relations: 
1) The average energy of a shwave grows in proportion to 
its frequency, and 2) the direction of angular momentum 
transport associated with a shwave can be determined 
from the time derivative of its frequency. The former 
implies that the energy of a compressive shwave (which 
has uj oc k or LO (X VAk when vax 0) grows asymptot- 
ically (i.e., when it is tightly wound) and is a minimum 
when it is in an open configuration (i.e., when kx = 0), 
whereas the energy of an incompressive shwave (which 
has UJ oc 1/fc) peaks at fca: = and decays asymptotically. 
The latter implies that a compressive shwave (which has 
to (X kx) transports angular momentum outward (inward) 
when it is in a trailing (leading) configuration, whereas 
an incompressive shwave (which has lu cx —kx) has the 
opposite behavior. ^2 



gives a comprehensive discussion of wave energy 
and angular momentum and their relation to action conservation 
in terms of spatial fluxes of definite frequency modes. There are 
no net spatial fluxes in a local model such as the shearing box, but 
one can determine the direction of global transport from the sign 
of the local shear stress, where a positiv e (negative) shear stress 
implies outward (inward) transp ort (e.g. JRvu &: GoodmanI 119921 : 
[Balbus & Hawlcv 1998; Balbus 2000h. 

These considerations assume that uj is real, a necessary condi- 
tion for the conservation of wave action; magnetorotationally un- 
stable modes are incompressive and yet transport angular momen- 
tum outwards. 
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6. SUMMARY AND DISCUSSION 

To within a multiphcative (dimensional) constant, the 
sohition for incompressive shwaves in an isothermal 
MHD shear flow is given to first order by: 



Sv — oj 



k X (k X X 



+i 



sb^ exp ^— z J ujdt + ik-x^, (80) 



5B 



Lj va ■ k 
. 2flkzUj 



k X (k X X 



exp ^— z J u! dt + ik ■ x^ , (81) 



and 



2n fk 



5p^-VA-5B~i— [ -^y+[q-l\^x]-5v, (82) 



k \ k 



where uj and uj are given by expressions ([28|) and ([29]) . 
The incompressive shwaves have 5v _L _L A;. They 
are not shear Alfvenic, since the effect of the shear is 
to increase the angle between va and k so that they 
are asymptotically perpendicular. Their energy peaks at 
fca; = and decreases as \kx\ oo, and they asymp- 
totically transport angular momentum inward. In gen- 
eral, these solutions are accurate for kz/ky 3> 1 (with 
errors < 30% for kz/ky < 1). Within the approxi- 
mations required for the assumption of incompressibil- 
ity {(3^1 and Cgk/n ^ 1), th ey are non linear so- 
lutions to the basic equations ( Goodman fc X u 1994). 
For wavevectors k initially aligned nearly perpendicu- 
lar to the mean magnetic field va, as well as for a very 
weak field {y/P :s> Cgk/il ^ 1), the solution correspond- 
ing to the lower branch of the dispersion relation ([^5]) is 
no longer uniformly valid. In that case, the solution is 
given by expressions ([50 |) -([56 |) multiplied by the factor 
exp {ik ■ x). 

To within a multiplicative (dimensional) constant, the 
solution for short-wavelength compressive shwaves 
is given to first order by: 



Sv 



5B 



ujk 

LO 



uj-'k 



j^-{vA-kY clB \k 



UJ 



k — Va ■ kvA 



X exp 



UJ dt + ik ■ X 



(83) 



uj^~{vA-k)clk^ 



VA-k 



X exp ^— i J UJ dt + ik ■ QC^ , (84) 



and 



J. ujk 



Cs V uj'^ - ivA-k)c^k^ 
X exp ( — i / UJ dt + ik ■ X 



(85) 



where here uj is given by expression (|62p . For the com- 
pressive shwaves, 5v and SB lie in the va, k plane, with 



Sv _L SB only for va -L k. Their energy is a minimum 
at fcj; = and increases as \kx\ — > oo, and they asymp- 
totically transport angular momentum outward. In gen- 
eral, these solutions are accurate for Hk ^ 1. The slow- 
shwave solution, corresponding to the lower branch of the 
dispersion relation (|62|) . breaks down for ^ 1, whereas 
the fast-shwave solution is generally valid for arbitrary 
p. Both solutions, however, can be compromised by the 
presence of a transiently-unstable incompressive shwave, 
and both break down when va \\ k. 

The solutions described here are ideal for convergence 
testing of MHD algorithms that incorporate a back- 
ground shear flow. They include both super- and sub- 
thermal magnetic fields as well as compressive and in- 
compressive velocity perturbations and have been veri- 
fied to be excellent approximations to the full set of linear 
equations for the bulk of wavenumbers that fit within the 
disk. Wave action is conserved for both sets of solutions, 
and is intimately related to the evolution of their energy 
and angular momentum transport (©. The outward (in- 
ward) transport of trailing compressive (incompressive) 
local disturbances is a general result. 

The methods employed in this paper can be used to 
extend any axisymmetric modal analysis in the shear- 
ing box to nonaxisymmetry. If the linear equations can 
be expressed in a form that does not depend explicitly 
on ky, the nonaxisymmetric dispersion relation is sim- 
ply given by the axisymmetric dispersion relation with 
the two-dimensional wavenumber replaced by the three- 
dimensional one, provided the resulting time-dependent 
frequency satisfies the adiabatic approximation. If, in ad- 
dition, the perturbation Lagrangian is known, the slowly- 
varying amplitude of the nonaxisymmetric modes can 
be calculated from the wave action of the axisymmetric 
modes. 

One possible astrophysical application is to study the 
nonlinear interaction of these s olutions or their stabi l- 
ity to secondary perturbations ([Goodman fc Xulll994f) : 
such a study might shed some light on the saturation 
of th e magnetorotational instability (Balbus fc HawlevI 
19911) in shearin g-box simulations (Hawlev et al. 1993 
Sano et al.ll2004[ ). In addition, the IGoldreich fc Sridhlirl 
(il995i) theory of MHD turbulence assumes that en- 
ergy injection at the forcing scale of the turbulence is 
roughly isotropic. Whereas the mode interactions in 
the Goldreich-Sridhar (GS) cascade take place on a time 
scale {vAk)~^ ^ and the bulk of the inertial range 
should therefore be unaffected by the shear, the incom- 
pressive modes considered here with VAk ^ fl have 
Va -L k for much of their evolution. It is not clear what 
effect energy injection into modes with va -L k will have 
upon the inertial range in the GS cascade (which consists 
of Alfven wave packets propagating along the field) , but 
it may be that measurements of the turbulent cascade 
in shea ring-box simu lations wi ll differ from the isotropic 
results (ICho fc Vishniac 200(31 : iMaron fc Goldreichll2001l: 
iGho et al.ll207il i! 
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APPENDIX 

A. HIGHER-ORDER EXPANSION FOR INCOMPRESSIVE SHWAVES 
The first-order equations for the incompressive shwaves are 

Sr](") - iuSr]'^'^^ - i{vA ■ k)k^Sv^Al + '^^k.Suji^^ = 0, (Al) 

SJ°^ - iujSoj^}^ -i{vA- k)64^^ - (2 - q)nk, ^ = 0, (A2) 
Svf, - iu5v^ll ^i{vA-k)^4i^^0, (A3) 



and 

S4°'> - tLjSji^'' - i(vA ■ k)Suj^^^ - qnkjv^2 = 0, (A4) 

where {3^°^} are the lowest-order solutions from ij3.2l and {(5^^^} are higher-order corrections. In the absence of the 
time derivatives of the zeroth-order amplitudes, the above set of equations is the same eigenvalue problem solved in 
il3.2l As a result, the combination of the first-order equations results in a term proportional to the dispersion relation 
(|28|) times one of the first-order amplitudes. This term is therefore of a higher order in the expansion and can be 
neglected, and one ends up with an equation for the time evolution of the zeroth-order amplitudes: 

■ (0) VA-k ■ (Q) mvA-kkz^- (0) nk^u ( oj^x . (o) 
OLOx oJx -I OVAx + ^ o 2 - df-j =0. (A5) 



Expressing all the amplitudes in terms of Sv'^ax eigenvector relationships (jSO]) and ((3T|) . which can also be 

written as 



and combining terms gives 



{2Cj'k' 



K k,) ^^[indvA,^j+ + 2^ ^^+2^ -0. (A7j 

Notice that this is a first-order differential equation; all of the degrees of freedom are contained within the dispersion 
relation ((28|) . 

Using the following relation obtained from differentiation of the dispersion relation ([28|: 

(fc^K^ „ 2(^2^2) ^^~4^ (AS) 

^ ' dt dt ^ ' 

equation (|A7p can be expressed as 

d / , (o)\ 1 d{Q^P) 1 ddi 1 duj ^ 

the direct integration of which yields expression (j37p . 

B. HIGHER-ORDER EXPANSION FOR COMPRESSIVE SHWAVES 
The first-order equations for the compressive shwaves are 

5p^°^ - iwJp(i) -f ik ■ Sv'^^^ = 0, (Bl) 

k ■ -t- ik^ {clSp'^^^ -f VA ■ (5?;^^^) + (2 - q)nkySv^°^ - 2nkjv^"'> = 0, (B2) 

VA ■ Sv^°'^ - iuJVA ■ Sv'-^^ + ivA ■ kclSp^^^ + (2 - q)VLvAy5v':^^ - 2nvAxSv^°^ = 0, (B3) 

and 

va ■ Sva — iujVA ■ Sv'"^ — ivA ■ kvA ■ Sv^^^ + iv\k ■ Sv'^^'' + qVlvAy^v'^Ax = 0, (B4) 

where {(5*^°^} are the lowest-order solutions from ij4.2l and {<5*-^^} are higher-order corrections. As in Appendix [Xl these 
can be combined into a single equation for the time dependence of the zeroth-order amplitudes: 

2/2^2 72 

k ■ 5v^^ -f Sp^"^ + {2- q)nkyx ■ 5v^"^ - 2nk^y ■ Sv^"^ + — (qilvAyX ■ 5v'^^'> + va ■ Sva''^^ 



k ■ ov — iLu 



VA ■ k 



VA ■ 5v'^^ + (2 - q)^VAyX ■ 5v^'^^ - 2nvAxy ■ 5v^^^ ) = 0. (B5) 



Notice that the time dependence of the mean magnetic field (equation [12]) imphes 

d 



— (vA ■ Sv'-"^^ = VA ■ 6v^°^ - qVlVAxV ■ 



(B6) 



with an analogous expression for va ■ ^^a ■ 

Considerable simplification of equation (|B5|) results from combining the radial and azimuthal components of the 
velocity and magnetic field, which are given by the dot product of x and y with ([55)1 and ([M]) . In particular, the 
Coriolis terms cancel, and the terms proportional to gil can be expressed in terms of derivatives of k and va ■ Performing 
this simplification and expressing all of the amplitudes in terms of k ■ (Jv^^^ via the eigenvector relationships and 
(l64l). rewritten as 



VA 



■5v^°\vA-5vf 



1, 



[va ■ k] 



one obtains the following equation for the amplitude of the compressive shwaves: 



(B7) 



C2fc2 {vA ■ kf d 



dt 



lnfe-5t;(°) 









1 UJ 


dt 


K ujk^ J 



{vA-kfcie d ( 1 







d 








+ 

W3 


dt 


\uj) 





u2 dv\ 
dt 



(B8) 



Differentiation of the dispersion relation with respect to t yields the following relationship between the time 
derivatives of fc, va and w: 

„2 , .2^U2 -[vA-kf duj-' 



dt B dt 



dt 



Using this and 



k d f bj 



w2 dt \ fc2 



fc2 d_ 

UJ^ dt \ 



[va ■ kf' c^k"^ 
k^ 



2 dio 
UJ dt ' 



and combining terms proportional to w, the amplitude equation (|B8p simplifies to 



d 
dt 



(ink-Sv'^^^ 



1 diu 3 duj 



[va ■ kf clk"^ 



Co dt 2uj dt uji-{vA-kfclk^dt 
Direct integration of the above expression gives 

fe • Sv'^^^ cx ujujk. 



(B9) 
(BIO) 

(Bll) 
(B12) 



Lo'^ — c-^k^ {va ■ k) 
which, along with k ■ Sv^'^'^ = uj6p^^\ is equivalent to expression ([66]) . 

C. SHWAVE ENERGY 

The second-order Lagrangian^'^ for a rotating, isothermal MHD shear flow (with po = 1) is given by pewaii 119701 : 
lHavea,1970; .Naravan et al..il987i) 



L2 = ^ 





di 




dt 



' dxj dxi 



dxj dxi 



nV-i] VA-[VA-Vi]-[VA-Vi\ 



(CI) 



Variation of this perturbation Lagrangian yields the linearized equations of motion (jl3p - (|16p . Inserting the real part 
of the perturbations (including the exponential factor) into expression (jCll) and averaging over azimuth gives 



2 [va • fe] [fe • i,VA ■C. + k- ^,VA ■ U - [va ■ fe] 



(C2) 



where |^|2 — ^'^ + S^f, and and ^^ are the real and imaginary parts of the fluid displacements (without the exponential 
factor). Noting that the fluid displacements are related to the density and magnetic field perturbations by 

Sp = -ik-£, (C3) 

Second-order here refers to an expansion in powers of a small-amplitude displacement. 
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and 

SB = i{vA-k^-k-ivA) (C4) 

via the Hnearized continuity and induction equations, the averaged perturbation Lagrangian can be expressed in simpler 
form as 

The average shwave energy as defined in expression (|34p is thus 

£^\{^' lel' - 2ql72 + c2 + \Sva\') . (C6) 
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Fig. 1. — Error in the energy of an incompressive shwave corresponding to the positive branch of the dispersion relation (|28^ for a series 
of numerical integrations with Cskz / (2-110.) = 10 and ky/k^ varying from 1 (top) to 10~* (bottom) by factors of 10. The data have been 
smoothed for readability. The shear parameter q = 3/2, corresponding to a Keplerian rotation profile. 




Fig. 2. — Error in the energy of an incompressive shwave corresponding to the negative branch of the dispersion relation (|28p for a series 
of numerical integrations with q = 3/2, Csky/{2iTQ) = 1 and kz/ky = 3 (top), 10, 30, 100 and 300 (bottom). The data have been smoothed 
for readability. 
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Fig. 3. — Error in the energy of a fast shwave for a series of numerical in tegr ations with q = 3/2 and Csky / {2itQ) = Cskz/{2TTQ) = 1 
(top), 3, 10, 30 and 100 (bottom). The scale factor / is defined in expression I I76II . The data have been smoothed for readability. 




Fig. 4. — Error in the energy of a slow shwave for a series of numerical in tegr ations with q = 3/2 and Csky / {2ttQ) = Csk^ / {2itVI) = 1 
(top), 3, 10, 30 and 100 (bottom). The scale factor / is defined in expression I I76II . The data have been smoothed for readability. 



